clc;
clear;

r = 5; % 阶数
n = 0.5*(r+1)*(r+2); % 节点数(基函数个数)

phi = cell(n,1);
node = zeros(n,2);
flag = 1;
for i = 0:r
    for j = 0:r
        if i+j <= r && i+j >= 0
            k = r-i-j;
            phi{flag,1} = func_N(i, j, k, r);
            node(flag,:) = [i, j]/r;
            flag = flag + 1;
        end
    end
end
disp(phi);

figure("WindowStyle", "docked");
title(["Lagrange Element Order r = " num2str(r)]);
scatter(node(:,1), node(:,2), "r*");
text(node(:,1)+0.005, node(:,2), int2str((1:n)'), "FontSize", 10, "FontWeight", "bold", "Color", "k");
axis("equal", "tight");

% 导数
syms x y lambda(x,y)
phi_dx = cell(n,1);
phi_dy = cell(n,1);
for i = 1:n
    phi_dx{i,1} = subs(diff(phi{i,1}, x), diff(lambda(x,y), x), -1);
    phi_dy{i,1} = subs(diff(phi{i,1}, y), diff(lambda(x,y), y), -1);
end
disp(phi_dx);
disp(phi_dy);



%% func_N
function N = func_N(i, j, k, r)
syms x y lambda(x,y)
t1 = sym(1);
for l = 0:i-1
    t1 = t1 * ((l-r*x)/(l-i));
end
t2 = sym(1);
for m = 0:j-1
    t2 = t2 * ((m-r*y)/(m-j));
end
t3 = sym(1);
for n = 0:k-1
    t3 = t3 * ((n-r*lambda(x,y))/(n-k));
end
N = t1*t2*t3;
end